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We demonstrate an unexpected connection between isotropic turbulence and wall-bounded shear 
flows. We perform direct numerical simulations of isotropic turbulence forced at large scales at 
moderate Reynolds numbers and observe sudden transitions from chaotic dynamics to a spatially 
simple flow, analogous to the laminar state in wall bounded shear flows. We find that the survival 
probabilities of turbulence are exponential and the typical lifetimes increase super-exponentially 
with the Reynolds number. Our results suggest that both isotropic turbulence and wall-bounded 
shear flows qualitatively share the same phase-space dynamics. 
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Recent years have seen significant advances in our 
understanding of the transition to turbulence in wall- 
bounded shear flows. In simple geometries, like flow in 
a pipe or a channel, close to the transition threshold, 
a finite-amplitude perturbation develops into a localised 
turbulent patch (a ‘puff’ in pipe flow) that exists as an 
independent entity m ■ Experiments pnni and nu¬ 
merical simulations nu-na have shown that the localised 
patches of turbulence can spontaneously disappear (re- 
laminarise) or split into two. The rates of these compet¬ 
ing processes depend strongly on the Reynolds number: 
at relatively low Reynolds numbers it is much more prob¬ 
able for a puff to decay than to split, while the opposite 
is true at higher Reynolds numbers. The point where 
the two probabilities are equal marks the transition to 
a sustained turbulence [TO], and the turbulence below 
this threshold may consist of long-lived chaotic transients 
m- The transition to turbulence in wall-bounded flows 
is thus intimately related to the process of relaminari¬ 
sation, where turbulent dynamics suddenly collapse to 
a much simpler, typically linearly stable, laminar state. 
Such events have been explained by dynamical systems 
theory as the escape from a chaotic saddle in state space 
with a constant (time independent) rate of escape [151- 
n. At higher Reynolds numbers, spatially local relam¬ 
inarisation attempts nano can be the source of inter- 
mittency in turbulent flows. 

In contrast, stationary isotropic turbulence, that can 
be thought of as a turbulent flow far away from bound¬ 
aries f22], is believed to exhibit much simpler dynam¬ 
ics: its motion is turbulent for all Reynolds numbers and 
there is no actual transition. In this Letter we report 
an unexpected connection between these two fields. We 
perform direct numerical simulations (DNS) of station¬ 
ary isotropic turbulence at low Reynolds numbers and 
observe sudden breakdowns of the turbulent dynamics 
in favour of a much simpler state. Similar observations 
have been made in connection to symmetry-breaking in 
isotropic turbulence [23] and in magnetohydrodynamic 
flows subject to electrical forcing [24] . We study the na¬ 
ture of this process and show that it is analogous to the 


relaminarisation events in wall-bounded parallel shear 
flows. We find that forced isotropic turbulence at rel¬ 
atively low Reynolds numbers is transient and the rate 
of its collapse is constant in time, resulting in exponen¬ 
tially distributed lifetimes of the turbulent state similar 
to pipe (9j [lDl [HI [17| and plane Couette flow puna in]. 

We perform direct numerical simulations of the incom¬ 
pressible Navier-Stokes equations (NSE) 

dtU = —VP — u • Vu + vAu + /, (1) 

V • u = 0, (2) 

where u denotes the velocity field, / is an external 
force, v is the kinematic viscosity, P is the pressure, 
and we set the density to unity. These equations were 
solved numerically using the standard fully de-aliased 
pseudospectral method m on a 3D periodic domain of 
length Lb 0 x = 27t with the smallest wavenumber being 
k m in = 27t/ L^ox = 1- All simulations are well-resolved, 
using 32 3 collocation points and satisfying k max ri ^ 1.82, 
where denotes the Kolmogorov dissipation scale. 

The system is forced at large scales by a negative 
damping f defined as 

f(k,t) = (s\y/2Ef)u(k,t) for 0 < |fc| < kf, 

= 0 otherwise. (3) 

Here, f(k,t) is the Fourier transform of the forcing, 
u(k , t ) is the Fourier transform of the velocity field 
u(x,t ), Ef is the total energy contained in the forcing 
band, and kf = 2.5 is the highest wavenumber forced. 
Normalizing the energy input by Ef ensures that the 
energy injection rate is £w = constant; here we choose 
£w = 0.1. This forcing provides an energy input that 
does not prefer any particular direction and has a compli¬ 
cated, time-dependent spatial profile; note that kf = 2.5 
corresponds to 80 possible wavevectors and thus 80 dif¬ 
ferent velocity field modes are being forced. It is com¬ 
monly used in numerical investigations of homogeneous 
isotropic turbulence [28ti34] , the prime example being the 
series of high-resolution simulations of Kaneda et al. ,35] • 
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The initial conditions for the velocity with a prescribed 
energy spectrum are constructed by assigning a Gaussian 
random vector to each point in space. The resulting field 
is subsequently Fourier-transformed and rescaled accord¬ 
ing to the desired energy spectrum in the form 

E(k) = 0.001702 fc 4 e" 2(fc/5)2 . (4) 

Further details of the numerical method, validation and 
benchmarking of the code can be found in m- 

The simulations are evolved for 1271 initial large-eddy 
turnover times to = L/U , where U denotes the initial 
rms velocity and L is the initial integral length scale; 
to = 0.78 in simulation units. The parameter that is var¬ 
ied in our simulations is the viscosity v , and we present 
the results in terms of a system-scale Reynolds number 
Re = L^e^ 3 /is that is changed from 53.80 to 97.82 for 
different simulations; in each individual run, Re is kept 
constant during the whole simulation. The corresponding 
values of the Taylor-Reynolds number and further simu¬ 
lation details are given in Supplementary Material [3B]. 

As mentioned above, the form of the forcing term that 
we employ here, Eq.([3]), is routinely used in DNS of 
isotropic turbulence as its complicated spatial form would 
seem to guarantee that the system is turbulent at any 
Reynolds number larger than unity. Indeed, even at suf¬ 
ficiently low Reynolds numbers, we observe that our sim¬ 
ulations reach a turbulent stationary state, where the en¬ 
ergy injection is balanced by the average dissipation and 
there is motion at all lengthscales. Surprisingly, however, 
after staying in this steady state for a long time, the sys¬ 
tem exhibits a transition to a different state, as shown, for 

example, in Fig.fTlfor Re = 76.86. There we plot the total 

1 — 1 £ 

energy of the system, E(t ) = J k max dk E(k , t ) and the en¬ 
ergy content of small scales, E'(t) = . dk E(k,t), 

as a function of time; the energy content at a particular 
lengthscale is E(k,t) = \u{k,t)\ 2 /2 and the largest scale 
in the system corresponds to k m i n = 2 tt/Li> ox = 1. 

As Fig. [l] demonstrates, turbulent dynamics persists 
until about t/to ~ 240. After that, the total energy be¬ 
comes constant and the small-scale fluctuations in the 
kinetic energy produced by the characteristic turbulent 
cascade process suddenly disappear. This implies that for 
t/to > 240 the kinetic energy is confined to the largest 
scale of the system and no nonlinear transfer exciting the 
smaller scales takes place. The system thus transitions 
from a turbulent to a large-scale ‘laminar’ state. 

The existence of such a state can be understood if one 
considers a model velocity field with u x ~ cos (y) and 
all other components of the velocity being zero. This 
flow profile is similar to a simple shear flow: it satisfies 
the incompressibility condition, it does not produce any 
pressure gradient in the system, and the non-linear term 
vanishes exactly for this profile. It is, therefore, an ex¬ 
act solution of the equations of motion, Eqs. 0 § , with 
its magnitude being set by the injection rate £w and 



FIG. 1. (Color online) Time evolution of the total energy E(t) 
and the energy content of the small scales E'(t) for Re = 76.86 
normalised by the initial energy Eq. Time is given in units 
of initial large eddy turnover time to = L/U, where U is the 
initial rms velocity and L the initial integral scale. The point 
around t/to ~ 240 when E'(t) vanishes and the total energy 
becomes constant marks the onset of the self-organised state 
as discussed in the main text. 

the kinematic viscosity is. In general, one can construct 
many exact solutions of the Navier-Stokes equations with 
k = 1, similar to the model profile discussed above, for 
which the non-linear term vanishes. What is surpris¬ 
ing, however, is that this self-organised large-scale state 
is dynamically connected to the isotropic turbulence at 
sufficiently low Reynolds numbers. 

When the system selects this self-organised state, it 
stays there for as long as our simulations continue. To¬ 
gether with the fact that this state is dynamically se¬ 
lected by the system, it seems to imply that this state 
is linearly stable. In order to further probe this state¬ 
ment, we have performed exploratory simulations where 
we have perturbed the self-organised state with random 
perturbations and observed their evolution. For suf¬ 
ficiently small amplitude of the perturbations, simula¬ 
tions always returned to the self-organised state, while 
for larger perturbations the system became turbulent, 
as shown in Fig. 2 of the Supplementary Material jMj . 
Therefore the simple state reported here has the same 
property as the laminar state in many wall-bounded par¬ 
allel shear flows (cf. the Hagen-Poiseuille profile in pipe 
flow |37]): it is a linearly stable simple exact solution that 
can be destabilised by a finite-amplitude perturbation. 

Next, we observe that at a fixed Reynolds number, 
the time of self-organisation ( t/to ~ 240 in the exam¬ 
ple above) strongly depends on the initial conditions. 
We explore this variability systematically by starting 
100 runs with different initial conditions for a fixed 
value of Re. In each simulation, we monitor the time- 
evolution of the total energy E(t) and the dissipation 
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FIG. 2. (Color online) Survival probability as a function of the 
dimensionless time i/to from the beginning of a simulation. 

rate e(t) = 2 u J^ max dk k 2 E(k,t). In order to identify 
the moment when the turbulent dynamics collapses onto 
the self-organised state, we employ a criterium that is 
based on the observation that since the kinetic energy in 
the self-organized state is confined to modes with k = 1, 
the asymptotic value E^ for all individual runs in a given 
ensemble (at a given Re) can be calculated from the en¬ 
ergy input rate £w and v. For statistically stationary 
flows the energy input rate Ew must equal the dissipa¬ 
tion rate e, and we obtain for the total energy of the 
self-organised state 

Eoc = £ -^ = constant. (5) 

Our data confirms that in every simulation, the total 
energy eventually reaches the asymptotic value Eoq , and 
the self-organisation time can be defined as the time when 
E(t) = £?«,. We have checked that when we define the 
relaminarisation moment as the time when the dissipa¬ 
tion rate equals the input rate without any fluctuations, 
we obtain identical results. 

We quantify the variability of the self-organisation 
times by introducing a survival probability P^ e (i) that 
at a given Re gives the probability that the system is still 
turbulent at time t, having started in a turbulent state 
at time t = 0. For each t, we estimate this probability 
by dividing the number of runs that are still turbulent 
after time t by the total number of runs performed at 
this Reynolds number. The resulting survival probabili¬ 
ties are shown in Fig(2]for a range of Re. We find that 
after some initial lag time during which the system has 
evolved from the initial condition into the turbulent state, 
the survival probability follows a simple exponential law 

Pne(t) ~ exp (-t/r(Re)), (6) 

where r(Re) is the typical lifetime of turbulence that 
only depends on the Reynolds number. The exponential 


form of the survival probability suggests that the process 
is memoryless, i.e. at each time the rate of relaminari¬ 
sation is constant and does not depend on the previous 
dynamics of the system. This behaviour is identical to 
what was observed in wall-bounded shear flows, such as 
pipe [HI m m HZIES] or plane Couette flow m mm- 
There, it was attributed to the escape from a chaotic 
saddle associated with relaminarisation of localised tur¬ 
bulence PUSHES]. 

In order to verify that our results do not depend on 
the size of the simulation box, one ensemble of 100 runs 
was created using a larger simulation box with Lb ox = 
47t. The collapse of turbulence is also observed in these 
runs and leads to an exponential survival probability with 
the same characteristic lifetime as a reference dataset at 
L box = 27T [SB]. 

The characteristic lifetime r is obtained at each 
Reynolds number from fitting the survival probabilities 
to Eq. ©>> see solid lines in Fig. [2j We observe a steep 
increase in r with increasing Reynolds number as shown 
in Fig. [3j To find the functional form r = r(Re), we fit 
the observed lifetime to various model expressions. First 
we consider a power law with an exponent n < 0 in the 
form r ~ (Re c — Re) n that would suggest a divergence of 
the lifetime at some critical Reynolds number Re c . We 
find that it is not compatible with the data for any value 
of n; Fig. 0 shows an example with n = —1. The same 
applies to an exponential increase of r with Re. However, 
we find that a super-exponential scaling in the form 

= cexp [exp(a + bRe)] (7) 

to 

is compatible with our data for a fixed amplitude c = 
15.63 and a = -3.48 ± 0.51, b = 0.052 ± 0.005, see 
Fig- 0 Once again, this conclusion parallels the super¬ 
exponential scaling of the lifetimes in wall-bounded shear 
flows (01 HH ESj • Further support for this scaling is pro¬ 
vided in the Supplementary Material [55] - 

We note that the super-exponential law Eq.([7]) is not 
the only possible form that produces an acceptable fit 
to the data. Another super-exponential dependence, 
r(Re)/to = exp[— a' + (6'-Re) 56 ]) with a' = —3.18 ± 0.14 
and b' = 0.0136 ± 0.0003 also gives a good agreement 
with the dataset, as can be seen in Fig. 0 

The results presented in this Letter show that there 
is a surprising analogy between the behaviour of the 
isotropic turbulence forced at a large scale and wall- 
bounded shear flows at low Reynolds numbers. We ob¬ 
serve that there is a spontaneous transition from turbu¬ 
lence to a spatially-simple state, which we have identified 
here, and this “laminar” state is linearly stable but can 
be destabilised by a finite-amplitude perturbation. The 
turbulent-laminar transition is abrupt and memory less, 
and the associated survival probability is exponential in 
time, cf. 0 HH EH HI ED EH- The turbulent lifetimes 
do not diverge with an increase in Re, but instead grow 




4 



Re 


FIG. 3. (Color online) Reynolds number dependence of the 
escape rate to/r. The red (grey) line is a two-parameter fit 
of the expression to/r(Re) = 0.064 exp (—exp [a + bRe]), the 
black line a two-parameter fit of the expression to/r(Re) = 
exp[a' — (b'Re) 5 b ]), the dash-dotted line a fit of an exponential 
and the faint dotted line a fit of a linear dependence of to/r 
on Reynolds number. 



Ei/Eo 

FIG. 4. Phase portrait E 2 vs E\ for Re = 76.86. Each point 
corresponds to a particular moment in time. All energies are 
scaled with the initial total kinetic energy Eq. Inset: Zoom 
of the turbulent region of the main graph showing that the 
dynamics is organised by several points in phase space sug¬ 
gestive of unstable exact solutions. 

super-exponentially, cf. [351 22- This analogy implies 
that the phenomena of the transition to turbulence in 
wall-bounded shear flows and forced isotropic turbulence, 
typically thought of as a high-Re phenomenon away from 
boundaries, are dynamically similar and can be under¬ 
stood within the same theoretical framework. As re¬ 
cent research suggests, the transition to turbulence in 
shear flows belongs to the directed percolation universal¬ 
ity class , and we argue that the same might 

be valid for forced isotropic turbulence. 

The phase space of turbulent wall-bounded shear flows 


is organised by exact solutions and periodic orbits of the 
Navier-Stokes equations [43 Hi! an d the relaminarisation 
events are associated with a sudden escape from this part 
of the phase space [TB]. Since we observe the same phe¬ 
nomenology, we speculate that the phase space of the 
forced isotropic turbulence should also be organised by 
coherent structures (exact solutions and periodic orbits). 
In Fig. [4] we plot the energy content in the k = 2 mode 
vs the energy in the k = 1 mode for a run at Re = 76.86. 
Each point there corresponds to a particular moment in 
time and the dynamics proceeds from left to right, until 
the system relaminarises (i.e. E\ = E 00 and E 2 = 0). We 
observe that the dynamics revolves around several points 
in phase space that are very suggestive of exact unstable 
solutions 51] ■ Identification of these coherent states will 
be the subject of future work. 

This work also suggests that the type of forcing em¬ 
ployed here is well-suited for DNS of isotropic turbulence 
with a view of creating an artificial, simpler system whose 
dynamics still resemble more complicated real physical 
systems, such as shear flows. Other types of forcing, no¬ 
tably various forms of stochastic forcing, are routinely 
used but might not have the phenomenological similar¬ 
ities with the transition to turbulence in wall-bounded 
shear flows. Recent results on self-organisation in mag¬ 
netohydrodynamic flows, for example, demonstrate that 
introducing random phases in the forcing term precludes 
the formation of a large-scale flow 123]. We argue that 
Eq.@ gives in fact a better approximation to naturally 
occurring turbulence than an explicitly stochastic (and 
thus more costly) forcing. 

Our results provide new potential targets for turbu¬ 
lence control, since we have shown that there is a stable 
large-scale state hidden in what appears to be isotropic 
turbulence. A particular choice of an additional external 
force may be sufficient to push the system into the basin 
of attraction of this stable state. 
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Supplemental Material 


Simulation parameters 

We perform direct numerical simulations of the incom¬ 
pressible Navier-Stokes equations, Eqs.(l) and (2) of the 
Main Text, with a forcing term using a fully de-aliased 
pseudospectral method m on a 3D periodic domain of 
length Lt ox = 2 tt. All simulations [JF are performed 
with 32 3 collocation points in a range of kinematic vis¬ 
cosities v changing from 0.1 to 0.055. This range corre¬ 
sponds to the Taylor-Reynolds number R\ = 2.61 — 4.72 
and the Reynolds number based on the box size (as de¬ 
fined in the Main Text) Re = 53.8 — 97.82; the prod¬ 
uct of the largest resolved wavenumber k max and the 
Kolmogorov lengthscale 77 is in the range of k max ri = 
2.85—1.82. Simulations are evolved for 1271t 0 time units, 
where the initial large-eddy turnover time to = L/JJ , 
with U and L being the initial rms velocity and integral 
lengthscale, correspondingly. For each Reynolds number 
we performed 100 runs starting from random initial con¬ 
ditions, as discussed in the Main Text. 


Stability of the self-organized state 

Here we test the stability properties of the self- 
organized state. Below, we consider an example of Re = 
75.78. Simulations are started from the self-organized 
state obtained at the end of a long run after a relam- 
inarization event at the same Reynolds number. This 
state is perturbed by a random initial perturbation of 
various amplitudes. In Fig|Sl|we plot the time evolution 
of three runs, with relatively small, medium and large 
amplitudes. We observe that the small and medium am¬ 
plitude runs return to the self-organized state after some 
transient dynamics, while the large amplitude run be¬ 
comes turbulent until its dynamics again collapses to the 
self-organized state. The lifetime of turbulence in the 
large-amplitude run is similar but not equal to the life¬ 
time of the original run (shown in Fig |Sl| for comparison), 
as can be expected from a process with an exponential 
survival probability (see Main Text). 


Survival Probabilities 

In order to further check that the statistics of relam- 
inarization events follow a simple exponential law with 
a super-exponential lifetime, we combine the results ob¬ 
tained in the Main Text for the survival probability yield¬ 
ing 



FIG. SI. Stability of the self-organized state for Re = 75.78. 


We now use the collection of the relaminarization times 
for various values of Re from our simulations to calculate 
the survival probability of the system being still turbulent 
after a fixed time t as a function of the Reynolds number. 
In Fig. |S2| we compare these results with the prediction 
of Eq.(Sl) for various values of the observation times t, 
where we held the constant dividing Re fixed while let¬ 
ting the additive constant a = 3.48 ± 0.51 vary within its 
error bounds calculated from the fitting procedure spec¬ 
ified in the Main Text. We observe a good agreement 
between the two data sets that provides further support 
to the exponential form of the survival probability with a 
superexponential lifetime as proposed in the Main Text. 
The characteristic S-shape of the curves is very similar 
to the results in wall-bounded flows 



f — 35 


20 exp [exp(—3.48 - 


Re \ 
19.23 ' 


P{t) = exp 


(SI) 


FIG. S2. Reynolds number dependence of the survival prob¬ 
abilities at different dimensionless observation times t/to- 
















Simulations using a larger box size 
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In order to demonstrate that the observed collapse of 
(transient) isotropic turbulence is not an artefact of a too 
small simulation box (L[, ox = 2n in the runs presented in 
the Main Text), an ensemble of 100 runs was created us¬ 
ing a box of size Lb ox = 47T. Figure [S3] shows the time evo¬ 
lution of the total energy E(t) and the energy content of 
the small scales E'(t ) for one simulation belonging to this 
ensemble. The collapse of the small-scale fluctuations is 
clearly visible, and the figure looks qualitatively very sim¬ 
ilar to Fig. 1 of the Letter. The survival probabilities for 
this ensemble and an ensemble using Lb ox = 2ir at the 
same Reynolds number are shown in Fig. |S4| where we 
note that the simulations using Lb ox = 4-7T take longer to 
reach a (transient) turbulent stationary state compared 
to the simulations using Lb ox = 27t. This is the reason for 
the shift between the survival probabilities visible in the 
figure. However, the slopes of the two exponentials, and, 
hence, the characteristic lifetimes obtained from the two 
datasets are almost indistinguishable, and certainly well 
within the error bars obtained from Eq. (SI). Since the 
simulations carried out on the larger box size show very 
similar features to the data obtained from simulation us¬ 
ing the conventional box size Lb ox = 27r, we conclude that 
the main results reported in the Letter are not induced 
by the size of the simulation box. 



r/r 0 



t/t 0 


FIG. S4. Survival probabilities as a function of the dimen¬ 
sionless time t/to from the beginning of a simulation using 
Lbox = 2-7T (black) and Lb ox = 47r (red) for Re = 68.10. The 
solid lines represent fits of P(t) = exp(—t/r) to the respective 
datasets. 


FIG. S3. Time evolution of the total energy E(t) and the 
energy content of the small scales E'(t) for Re = 68.10 nor¬ 
malized by the initial energy Eo for Lb ox = 47r. Time is given 
in units of initial large eddy turnover time to = L/U, where 
U is the initial rms velocity and L the initial integral scale. 














